function dvf_pseudo = create_psuedo_dvf_vec_volt3d_sliding_IterateVoltage_parfor(patient, volt, flow, scan, ii,jj,kk)
  %------------------------------------------------------------------------
    %   This file is part of the
    %   5D-Novel4DCT Toolbox  ("Novel4DCT-Toolbox")
    %   DH Thomas, Ph.D
    %   University of California, Los Angeles
    %   Contact: mailto:dhthomas@mednet.ucla.edu
    %------------------------------------------------------------------------
    % $Author: DHThomas $	$Date: 14-Apr-2014 15-23-24 $	$Revision: 0.1 $


tic
% mkdir('Psuedo_dvfs')

    dvf_pseudo=zeros(size(ii,2),size(jj,2),3,size(kk,2));
  tic 
%   cd(model_folder)
  
      volt_2d = permute(repmat(volt, [1 patient.dim(2)]),[2,1]);
    flow_2d = permute(repmat(flow, [1 patient.dim(2)]),[2,1]);
    k_2d = permute(repmat([1:patient.dim(3)]', [1 patient.dim(2)]),[2,1]);
    j_2d = permute(repmat([1:patient.dim(2)]', [1 patient.dim(3)]),[2,1])';

%     p_volt = polyfit(1:length(volt),volt',15);
%     p_flow = polyfit(1:length(flow),flow',15);
%     volt_1000 = spline([1:length(volt)],volt,[1:0.01:length(volt)]);
%     volt_1000 = polyval(p_volt, [1:0.0001:length(volt)]);
%     flow_1000 = polyval(p_flow, [1:0.0001:length(volt)]);
dvf_pseudo_u = zeros(patient.dim);
dvf_pseudo_v = zeros(patient.dim);
dvf_pseudo_w = zeros(patient.dim);

if patient.par_toolbox>0
    pctRunOnAll javaaddpath E:\ParforProgMonv2\java
    progressStepSize = 1;
    ppm = ParforProgMon('Creating Psuedo Deformation Vector: ',patient.dim(1),  progressStepSize,800, 120);
    
    parfor slice=1:size(ii,2);
        
        
        constant=load([patient.model_params_folder, sprintf('/constant_%d',slice)]);
        constant_vec_u=constant.data(:,:,1);constant_vec_u=constant_vec_u(:);
        constant_vec_v=constant.data(:,:,2);constant_vec_v=constant_vec_v(:);
        constant_vec_w=constant.data(:,:,3);constant_vec_w=constant_vec_w(:);
        
        alpha=load([patient.model_params_folder,sprintf('/alpha_%d',slice)]);
        alpha_vec_u=alpha.data(:,:,1);alpha_vec_u=alpha_vec_u(:);
        alpha_vec_v=alpha.data(:,:,2);alpha_vec_v=alpha_vec_v(:);
        alpha_vec_w=alpha.data(:,:,3);alpha_vec_w=alpha_vec_w(:);
        
        beta=load([patient.model_params_folder,sprintf('/beta_%d',slice)]);
        beta_vec_u=beta.data(:,:,1);beta_vec_u=beta_vec_u(:);
        beta_vec_v=beta.data(:,:,2);beta_vec_v=beta_vec_v(:);
        beta_vec_w=beta.data(:,:,3);beta_vec_w=beta_vec_w(:);
        
        %     gamma=load(sprintf('gamma_%d',slice));
        %     gamma_vec_u=alpha.data(:,:,1);gamma_vec_u=gamma_vec_u(:);
        %     gamma_vec_v=alpha.data(:,:,2);gamma_vec_v=gamma_vec_v(:);
        %     gamma_vec_w=alpha.data(:,:,3);gamma_vec_w=gamma_vec_w(:);
        
        model_u=[constant_vec_u alpha_vec_u beta_vec_u];% gamma_vec_u];
        model_v=[constant_vec_v alpha_vec_v beta_vec_v];% gamma_vec_v];
        model_w=[constant_vec_w alpha_vec_w beta_vec_w];% gamma_vec_w];
        
        E=[ones(length(model_u),1)  volt_2d(:) flow_2d(:)];
        dvf_fit_u=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
        dvf_fit_v=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
        dvf_fit_w=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
        
        %%
        
%         load([patient.dvf_folder sprintf('/dvfs_slice_%d',slice)]);
%         dvf_elastix_w = squeeze(dvf_tmp(:,:,:,3)); clear dvf_tmp
%         dvf_elastix_w_tmp = dvf_elastix_w(:,:);
%         %
%         [~,grid_z] = meshgrid(jj,kk);
%         grid_z=permute(grid_z,[2,1]);
%         voxel_w=double(bsxfun(@plus,grid_z,dvf_elastix_w));
%         voxel_w = voxel_w(:);
%         volt_2s = nakeinterp1([1:patient.dim(3)]',volt,voxel_w(:));
%         volt_2s = reshape(volt_2s,size(volt_2d));
%         flow_2s =  nakeinterp1([1:patient.dim(3)]',flow,voxel_w(:));
%         flow_2s = reshape(flow_2s,size(volt_2d));
%             
%         E_s=[ones(length(model_u),1)  volt_2s(:) flow_2s(:)];

            %      X3(scan,:) = nakeinterp1(kk',ecg_volt(kk,scan),voxel_w(scan,:)');
       
            %%
            % The bellows voltage in the 
        
        E_iterate = E;
 p_max =3;
 E_diff = zeros([p_max, size(E)]);
        for vox = 1:length(model_u); 
            %         if mod(vox,1000)==1
            %             vox
            %         end
%             vox_iterate = vox;
            for p = 1:p_max;
                
                dvf_fit_w(vox)=model_w(vox,:)*E_iterate(vox,:)';
                %             [j,k] = ind2sub(size(volt_2d),vox);
                %             volt_1d = polyval(p_volt, k_2d(vox)+(dvf_fit_w(vox)));
                %             flow_1d = polyval(p_flow, k_2d(vox)+(dvf_fit_w(vox)));
                if round((k_2d(vox)+(dvf_fit_w(vox))))>0&& round((k_2d(vox)+(dvf_fit_w(vox))))<length(flow);
                    volt_1d = volt(round((k_2d(vox)+(dvf_fit_w(vox)))));
                    flow_1d  =  flow(round((k_2d(vox)+(dvf_fit_w(vox)))));
                     
                    E_diff(p,vox,:) = E_iterate(vox,:)-[1 volt_1d flow_1d];
                    E_iterate(vox,:) = [1 volt_1d flow_1d];
                    
%                     vox_iterate = sub2ind([size(j_2d)'], j_2d(vox),round(k_2d(vox)+dvf_fit_w(vox)));

                end
            end
        end

        
        %     for vox = 1:length(model_u);
        %         dvf_fit_u1(vox)=model_u(vox,:)*E(vox,:)';
        %         dvf_fit_v1(vox)=model_v(vox,:)*E(vox,:)';
        %         dvf_fit_w1(vox)=model_w(vox,:)*E(vox,:)';
        %
        %     end
        for vox = 1:length(model_u);
            dvf_fit_u(vox)=model_u(vox,:)*E_iterate(vox,:)';%E_s(vox,:)';%
            dvf_fit_v(vox)=model_v(vox,:)*E_iterate(vox,:)';%E_s(vox,:)';%
            dvf_fit_w(vox)=model_w(vox,:)*E_iterate(vox,:)';%E_s(vox,:)';%
            
        end
        
        %     cd('Psuedo_dvfs');
        
        %     dvf_pseudo_u1(slice,jj,kk)=reshape(dvf_fit_u1,size(jj,2),size(kk,2));
        %     dvf_pseudo_v1(slice,jj,kk)=reshape(dvf_fit_v1,size(jj,2),size(kk,2));
        %     dvf_pseudo_w1(slice,jj,kk)=reshape(dvf_fit_w1,size(jj,2),size(kk,2));
        
        dvf_pseudo_u(slice,jj,kk)=reshape(dvf_fit_u,size(jj,2),size(kk,2));
        dvf_pseudo_v(slice,jj,kk)=reshape(dvf_fit_v,size(jj,2),size(kk,2));
        dvf_pseudo_w(slice,jj,kk)=reshape(dvf_fit_w,size(jj,2),size(kk,2));
        
        
        %          if save_files==1;
        %
        % %         iSaveX( sprintf( 'dvf_pseudo_%d', slice ), dvf_pseudo );
        %     end
        ppm.increment();
        
    end
else

    for slice=128;%:size(ii,2);
        display(sprintf('Scan %d; Generating DVF: Slice %d of %d', scan, slice, size(ii,2)));
        
        constant=load([patient.model_params_folder, sprintf('/constant_%d',slice)]);
        constant_vec_u=constant.data(:,:,1);constant_vec_u=constant_vec_u(:);
        constant_vec_v=constant.data(:,:,2);constant_vec_v=constant_vec_v(:);
        constant_vec_w=constant.data(:,:,3);constant_vec_w=constant_vec_w(:);
        
        alpha=load([patient.model_params_folder,sprintf('/alpha_%d',slice)]);
        alpha_vec_u=alpha.data(:,:,1);alpha_vec_u=alpha_vec_u(:);
        alpha_vec_v=alpha.data(:,:,2);alpha_vec_v=alpha_vec_v(:);
        alpha_vec_w=alpha.data(:,:,3);alpha_vec_w=alpha_vec_w(:);
        
        beta=load([patient.model_params_folder,sprintf('/beta_%d',slice)]);
        beta_vec_u=beta.data(:,:,1);beta_vec_u=beta_vec_u(:);
        beta_vec_v=beta.data(:,:,2);beta_vec_v=beta_vec_v(:);
        beta_vec_w=beta.data(:,:,3);beta_vec_w=beta_vec_w(:);
        
        %     gamma=load(sprintf('gamma_%d',slice));
        %     gamma_vec_u=alpha.data(:,:,1);gamma_vec_u=gamma_vec_u(:);
        %     gamma_vec_v=alpha.data(:,:,2);gamma_vec_v=gamma_vec_v(:);
        %     gamma_vec_w=alpha.data(:,:,3);gamma_vec_w=gamma_vec_w(:);
        
        model_u=[constant_vec_u alpha_vec_u beta_vec_u];% gamma_vec_u];
        model_v=[constant_vec_v alpha_vec_v beta_vec_v];% gamma_vec_v];
        model_w=[constant_vec_w alpha_vec_w beta_vec_w];% gamma_vec_w];
        
        E=[ones(length(model_u),1)  volt_2d(:) flow_2d(:)];
        dvf_fit_u=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
        dvf_fit_v=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
        dvf_fit_w=zeros(size(jj,2)*size(kk,2),1);%,size(ii,2),
        
        
               E_iterate = E;
 p_max =3; % perform 3 iterations to find the correct voltage.
 E_diff = zeros([p_max, size(E)]);


        for vox = 1:length(model_u);
            %         if mod(vox,1000)==1
            %             vox
            %         end
            vox_iterate = vox;
             for p = 1:p_max;
                
                dvf_fit_w(vox)=model_w(vox,:)*E_iterate(vox,:)';
                %             [j,k] = ind2sub(size(volt_2d),vox);
                %             volt_1d = polyval(p_volt, k_2d(vox)+(dvf_fit_w(vox)));
                %             flow_1d = polyval(p_flow, k_2d(vox)+(dvf_fit_w(vox)));
                if round((k_2d(vox)+(dvf_fit_w(vox))))>0&& round((k_2d(vox)+(dvf_fit_w(vox))))<length(flow);
                    volt_1d = volt(round((k_2d(vox)+(dvf_fit_w(vox)))));
                    flow_1d  =  flow(round((k_2d(vox)+(dvf_fit_w(vox)))));
                     
%                     E_diff(p,vox,:) = E_iterate(vox,:)-[1 volt_1d flow_1d];
                    E_iterate(vox,:) = [1 volt_1d flow_1d];
                    
%                     vox_iterate = sub2ind([size(j_2d)'], j_2d(vox),round(k_2d(vox)+dvf_fit_w(vox)));

                end
            end
        end
        
        %     for vox = 1:length(model_u);
        %         dvf_fit_u1(vox)=model_u(vox,:)*E(vox,:)';
        %         dvf_fit_v1(vox)=model_v(vox,:)*E(vox,:)';
        %         dvf_fit_w1(vox)=model_w(vox,:)*E(vox,:)';
        %
        %     end
        for vox = 1:length(model_u);
            dvf_fit_u(vox)=model_u(vox,:)*E_iterate(vox,:)';
            dvf_fit_v(vox)=model_v(vox,:)*E_iterate(vox,:)';
            dvf_fit_w(vox)=model_w(vox,:)*E_iterate(vox,:)';
            
        end
        
        %     cd('Psuedo_dvfs');
        
        %     dvf_pseudo_u1(slice,jj,kk)=reshape(dvf_fit_u1,size(jj,2),size(kk,2));
        %     dvf_pseudo_v1(slice,jj,kk)=reshape(dvf_fit_v1,size(jj,2),size(kk,2));
        %     dvf_pseudo_w1(slice,jj,kk)=reshape(dvf_fit_w1,size(jj,2),size(kk,2));
        
        dvf_pseudo_u(slice,jj,kk)=reshape(dvf_fit_u,size(jj,2),size(kk,2));
        dvf_pseudo_v(slice,jj,kk)=reshape(dvf_fit_v,size(jj,2),size(kk,2));
        dvf_pseudo_w(slice,jj,kk)=reshape(dvf_fit_w,size(jj,2),size(kk,2));
        
        
        %          if save_files==1;
        %
        % %         iSaveX( sprintf( 'dvf_pseudo_%d', slice ), dvf_pseudo );
        %     end
        %     ppm.increment();
        
    end
end

toc
dvf_pseudo(:,:,1,:)=dvf_pseudo_u;
dvf_pseudo(:,:,2,:)=dvf_pseudo_v;
dvf_pseudo(:,:,3,:)=dvf_pseudo_w;






